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Abstract 

This  report  describes  the  formulation  of  the  SEDZLJ  sediment  transport 
model,  coupling  of  the  model  with  the  ICM  eutrophication  model, 
validation  of  the  combined  codes  on  a  simplified  grid,  and  application  of 
the  combined  codes  to  a  prototype  system. 

The  code  validation  tests  indicate  the  codes  are  correctly  coupled.  The 
combined  codes  meet  specified  performance  criteria  including: 

1.  Mass  conservation  in  water  and  bed. 

2.  Limited  sensitivity  to  variations  in  model  time-step. 

3.  Settling  of  new  ICM  state  variables  through  model  water  column  agrees 
with  settling  of  original  ICM  state  variables. 

4.  Transport  of  new  ICM  state  variables  agrees  with  transport  of  original  ICM 
state  variables. 

5.  Sediment  bed  will  armor. 

6.  Bed  erodes  away  smoothly. 

7.  Model  behaves  reasonably  for  accumulation  in  the  bed  with  no  erosion. 

Additional  evaluation  and  testing  were  performed  through  application  of 
the  combined  models  to  a  prototype  system:  Lake  George,  Florida.  The 
codes  perform  reasonably  through  a  three-year  simulation.  The  tests 
indicate,  however,  that  realistic  application  of  SEDZLJ  requires  an 
extensive  set  of  field  and  laboratory  observations. 
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1  Introduction 

This  report  describes  the  merger  of  the  SEDZLJ  sediment  transport  model 
with  the  U.S.  Army  Corps  of  Engineers’  CE-QUAL-ICM  (or  simply  ICM) 
water  quality  model.  The  merger  of  ICM  with  SEDZLJ  is  the  first  step  in  the 
development  of  a  high-fidelity  contaminant  transport  and  fate  model  based 
on  the  ICM  concept  and  formulations.  The  sediment  transport  model  is 
necessary  because  contaminants  commonly  partition  to  sediments  of 
various  compositions  and  are  therefore  transported  along  with  the 
sediments.  Burial  in  bottom  sediments  provides  one  pathway  by  which 
contaminants  leave  a  system.  Quantification  of  this  removal  requires 
knowledge  of  long-term  sediment  erosion  and  deposition  rates.  The 
coupling  of  the  sediment  transport  model  with  ICM  takes  advantage  of  the 
detailed  ICM  representation  of  the  aquatic  carbon  cycle  in  the  water  column 
and  bed  sediments,  which  is  necessary  since  hydrophobic  contaminants 
display  a  strong  tendency  to  adhere  to  organic  carbon  particles. 

Numerous  sediment  transport  models  exist  (Warner  et  al.  2008,  Jones  and 
Lick  2001,  Tetra  Tech  Inc.  2007)  and  there  was  no  need  to  develop  a  new 
one  for  this  application.  SEDZLJ  was  selected  because  of  its  previous  use  in 
contaminant  fate  and  transport  studies  (e.g.  Weston  Solutions  Inc.  2006), 
because  Corps  personnel  have  experience  with  the  code  and  formulations, 
and  because  the  model  is  readily  parameterized  using  results  from  the 
Corps’  SEDFLUME  device. 

Anchor  QEA  LLC,  Montvale,  New  Jersey,  was  contracted  to  merge  the 
SEDZLJ  code  into  ICM.  Personnel  at  the  U.S.  Army  Engineer  Research 
and  Development  Center  (ERDC)  performed  quality  assurance  checks  on 
the  merger  and  completed  the  coupling  of  SEDZLJ  with  various  facets  of 
ICM.  This  report  describes  the  coupling  of  the  two  models,  basic 
performance  testing,  and  implementation  of  the  merged  models  on  a 
prototype  system. 
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2  Coupling  SEDZLJ  to  ICM 

State  variables 

Three  new  state  variables  were  introduced  to  the  ICM  suite: 


•  Clay/Silt  (Cohesive  Sediment) 

•  Fine  Sand  (Non-cohesive  Sediment) 

•  Medium  Sand  (Non-cohesive  Sediment) 


In  addition,  two  existing  ICM  variables,  labile  particulate  organic  carbon 
(LPOC)  and  refractory  particulate  organic  carbon  (RPOC)  were  treated  as 
discrete  cohesive  sediment  particles.  This  treatment  avoids  the  need  to 
assign  carbon  fractions  to  mineral  particles  as  in  classic  models  of 
contaminant  partitioning  (Karickhoff  et  al.  1979).  The  incorporation  of 
discrete  carbon  particles  allows  for  the  use  of  the  full  ICM  carbon  kinetics 
algorithms  and  allows  for  differential  transport  of  carbon  and  mineral 
particles. 


Conservation  of  mass  equation 

Foundation 


The  foundation  of  CE-QUAL-ICM  is  the  solution  to  the  three-dimensional 
mass-conservation  equation  for  a  control  volume.  Control  volumes 
correspond  to  cells  on  the  model  grid.  CE-QUAL-ICM  solves,  for  each 
volume  and  for  each  state  variable,  the  equation: 


g  Vj  ■  Cj 

8  t 


=  D?*  '  + 


Dk 


k=i 


k=i 


8  C 
8  xk 


+  2  Sj 


(1) 


where 

Vj  =  volume  of  jth  control  volume  (m3) 

Cj  =  concentration  in  jth  control  volume  (g  nr3) 
t,  x  =  temporal  and  spatial  coordinates 
n  =  number  of  flow  faces  attached  to  jth  control  volume 
Qk  =  volumetric  flow  across  flow  face  k  of  jth  control  volume  (m3  s1) 
Ck  =  concentration  in  flow  across  face  k  (g  m  3) 
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Ak  =  area  of  flow  face  k  (m2) 

Dk  =  diffusion  coefficient  at  flow  face  k  (m2  s1) 

Sj  =  external  loads  and  internal  sources  and  sinks  in  jth  control 
volume  (g  S'1) 

For  the  sediment  particles  in  the  water  column,  the  only  internal 
source/sink  is  settling: 


o 

S:  =  -Ws—SSI 
J  8z 


(2) 


where 

Ws  =  settling  velocity  (m  d  l) 

SSI  =  solids  concentration  (g  m  3) 
z  =  vertical  coordinate  (m) 

Settling  velocity  is  computed  by  Cheng’s  formula  (Cheng  1997): 


Ws 


v 

d 


■y/25+1.2  dl 


(3) 


where 


v  =  kinematic  fluid  viscosity  (cm2  s_1) 
d  =  particle  diameter  (cm) 

Parameter  d  is  non-dimensional  particle  diameter 


l  n2  J 


(4) 


where 


g  =  gravitational  acceleration  (cm  s  2) 


and 


_  Ps-  P 


A 


P 


(5) 
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where 


ps  =  density  of  particle  (g  cm  3) 
p  =  density  of  water  (g  cm  3) 


Settling  velocity  of  carbon  particles  is  a  model  input  variable. 

Deposition  and  erosion 

The  sediment  source/sink  term  for  ICM  cells  that  intersect  the  bottom  is 
modified  to  include  deposition  to  and  erosion  from  the  bed: 


Ws  E-D 


(6) 


where 

SSIup  =  solids  concentration  in  cell  overlying  cell  j  (g  cm -3) 
E  =  erosion  rate  (g  cm2  s1) 

D  =  deposition  rate  (g  cm2  s1) 

Az  =  cell  thickness  (m) 

The  erosion  rate  is  determined  by: 


E  =  fk-  BDD  ■  Sk  ■  Psusk  ■  Eg  r  oss 


(7) 


where 


BDD 

Sk 

Psusk 

Egross 


fraction  of  sediment  class  k  in  active  surface  layer  (o  <  fk  <  l) 
bed  dry  density  (g  cm  3) 

particle  shielding  factor  for  sediment  class  k  (o  <  Sk  <  l) 
probability  of  suspension  for  sediment  class  k  (o  <  Psusk  <  l) 
gross  erosion  rate  (cm  s1) 


The  gross  erosion  rate  is  computed  as: 

Egross  =  A- Tnsf  for  rsf  >  t„. 

=  0  f0rTsf  <  Tcr 


(8) 
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where 

TSf  =  bottom  shear  stress  (dyne  cm'2) 

Tcr  =  critical  shear  stress  for  erosion  (dyne  cm-2) 

A,n  =  empirical  parameters 

Particle  shielding  and  probability  of  suspension  are  described  in  Jones  and 
Lick  (2001). 

Deposition  is  computed: 


D  =  PdepkWskSk 


(9) 


where 

Pdepk  =  probability  of  deposition  (o  <  Pdepk  <  1) 

Probability  of  deposition  is  determined  for  cohesive  and  non-cohesive 
sediment  classes  as  described  by  Jones  and  Lick  (2001). 
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3  Basic  Performance  Tests 

Basic  performance  criteria  were  specified  for  acceptance  of  the  merger  of 
the  SEDZLJ  and  ICM  codes.  The  criteria  were  developed  based  on  initial 
explorations  of  the  linked  models  and  included: 

1.  Mass  conservation  in  water  and  bed. 

2.  Limited  sensitivity  to  variations  in  model  time-step. 

3.  Settling  of  new  ICM  state  variables  through  model  water  column  agrees 
with  settling  of  original  ICM  state  variables. 

4.  Transport  of  new  ICM  state  variables  agrees  with  transport  of  original  ICM 
state  variables. 

5.  Sediment  bed  will  armor. 

6.  Bed  erodes  away  smoothly. 

7.  Model  behaves  reasonably  for  accumulation  in  the  bed  with  no  erosion. 

The  explorations  were  conducted  on  a  30-box  grid  (Figure  1)  developed  as 
an  ICM  test  bed.  Geometry  and  circulation  in  the  test  system  are  roughly 
scaled  to  resemble  the  mainstem  of  Chesapeake  Bay.  The  test  bed  provides 
the  developer  with  maximum  flexibility  to  examine  model  behavior  via 
modifications  to  the  ICM  inputs  and  options  installed  in  the  ICM  code. 
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Figure  1.  Elevation  view  of  30-box  test  grid.  Cells  are  30  km  (length)  x  20  km  (width)  x  5  m 

(height). 


Mass  conservation 

For  this  test,  horizontal  transport  was  eliminated.  Mass  conservation  was 
examined  in  columns  consisting  of  three  water  cells  and  five  sediment 
layers.  Active  processes  included  vertical  diffusion,  settling,  erosion,  and 
deposition.  The  time-step  was  1000  s  and  the  test  duration  was  five  days. 


ER  DC/EL  TR-12-17 


7 


Initial  conditions  were  identical  for  all  constituents  and  provided  a 
uniform  distribution  in  the  water  column. 

Over  the  five-day  period,  the  carbon  components  exhibited  minimal 
settling  through  the  water  column  (Figure  2).  Concentration  increased  in 
the  bottom  water  layer,  however,  due  to  erosion  from  the  sediments.  The 
clay/ silt  component  exhibited  noticeable  settling  through  the  water 
column,  due  to  higher  settling  velocity,  and  increased  concentration  in  the 
bottom  water  layer  due  to  material  settling  from  above  and  eroding  from 
below  (Figure  3).  Fine  sand  settled  completely  out  of  the  upper  two  water 
layers  and  collected  in  the  bottom  layer  due  to  settling  from  above  and 
erosion  from  below  (Figure  4).  Medium  sand  disappeared  completely  from 
the  water  and  accumulated  in  the  sediments  (Figure  5).  Roughly  0.06%  of 
the  initial  mass  of  each  component  was  lost  over  the  five-day  period 
(432  time  iterations,  Table  1). 

Model  time-step 

The  mass  conservation  tests  were  repeated  with  two  different  time-steps: 
250  s  and  2000  s.  The  model  demonstrated  no  degenerate  behavior  or 
extreme  sensitivity  to  the  time-step  (Table  1).  Results  were  slightly  different 
for  each  test  with  mass  potentially  lost  or  gained.  Results  suggested  a  loss  of 
accuracy  as  the  time-step  was  increased  with  the  least  error  («  0.04%)  at  the 
250-s  step  and  the  greatest  error  (~  0.09%)  at  the  2000-s  step.  Accuracy 
was  slightly  less  for  medium  sand,  which  had  the  greatest  settling  velocity, 
than  for  POC  and  clay/silt,  which  had  the  slowest  settling  velocities.  The 
small  mass  gains  or  losses  were  attributed  to  numerical  round-off  error  and 
were  deemed  acceptable. 

Settling  velocity 

ICM  treats  particle  settling  as  a  term  in  the  kinetics  formulations  and  codes. 
To  ensure  that  the  settling  of  the  new  and  modified  constituents  agreed  with 
results  from  previously  tested  and  validated  code,  the  computed  concentra¬ 
tions  of  clay/silt  (new  constituent)  and  LPOC  (modified  constituent)  were 
compared  to  fixed  solids  (original  ICM  variable)  under  conditions  of  pure 
settling.  This  test  used  the  same  conditions  as  the  mass  balance  test  except 
erosion  and  deposition  were  eliminated.  Settling  velocity  of  the  three 
particle  types  was  set  to  1.319  m  d1  and  the  run  was  executed  for  10  days.  At 
the  end  of  the  run,  concentrations  of  fixed  solids,  LPOC,  and  clay/ silt  were 
identical  (Figure  6). 
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Figure  2.  LPOC  results  from  mass  balance  test.  Initial  and  final  concentrations  are  shown  for 
three  levels  in  the  water  column  and  for  the  surface  bed  layer.  Results  for  RPOC  are  identical. 
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Figure  3.  Clay/silt  results  from  mass  balance  test.  Initial  and  final  concentrations  are  shown 
for  three  levels  in  the  water  column  and  for  the  surface  bed  layer. 
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Figure  4.  Fine  sand  results  from  mass  balance  test.  Initial  and  final  concentrations  are  shown 
for  three  levels  in  the  water  column  and  for  the  surface  bed  layer. 
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Figure  5.  Medium  sand  results  from  mass  balance  test.  Initial  and  final  concentrations  are 
shown  for  three  levels  in  the  water  column  and  for  the  surface  bed  layer. 
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Table  1.  Mass  conservation  and  sensitivity  to  time-step  test  results.  Concentrations  (g  nr2)  are 
presented  for  a  single  column  of  water  and  sediment  cells. 


Variable 

Initial 

Concentration 

Final 

Concentration, 
At  =  250  s 

Final 

Concentration, 
At  =  1000  s 

Final 

Concentration, 
At  =  2000  s 

LPOC,  water 

150 

2,857.45 

2,870.80 

2,861.50 

LPOC,  sediments 

13,200 

10,495.94 

10,471.40 

10,479.35 

LPOC,  total 

31,350 

13,353.39 

13,342.20 

13,340.90 

RPOC,  water 

150 

2,857.45 

2,870.80 

2,861.50 

RPOC,  sediments 

13,200 

10,495.94 

10,471.40 

10,479.35 

RPOC,  total 

13,350 

13,353.39 

13,342.20 

13,340.90 

Clay/Silt,  water 

150 

2,857.45 

2,870.85 

2,861.55 

Clay/Silt,  sediments 

13,200 

10,495.94 

10,471.40 

10.479.35 

Clay/Silt,  total 

13,350 

13,353.39 

13,342.25 

13,340.90 

Fine  Sand,  water 

150 

2,907.55 

2,906.75 

2,826.40 

Fine  Sand,  sediments 

13,200 

10.446.17 

10,438.32 

10,518.06 

Fine  Sand,  total 

13,350 

13,353.72 

13,345.07 

13,344.46 

Medium  Sand,  water 

150 

0.00 

0.00 

0.00 

Medium  Sand,  sediments 

13,200 

13,354.95 

13,341.80 

13,338.36 

Medium  Sand,  total 

13,350 

13,354.95 

13,341.80 

13,338.36 

Transport 

One  of  the  goals  of  this  research  was  to  ensure  that  transport  of  new  and 
modified  variables  agreed  with  previously  tested  and  validated  code.  This 
test  was  similar  to  the  test  of  settling  except  that  two-dimensional  circula¬ 
tion  similar  to  long-term  average  circulation  in  Chesapeake  Bay  was 
imposed.  This  run  incorporated  longitudinal  and  vertical  currents,  longi¬ 
tudinal  and  vertical  diffusion,  and  vertical  settling.  At  the  end  of  a  10-day 
model  run,  computed  concentrations  of  fixed  solids,  LPOC,  and  clay/ silt 
were  virtually  identical  (Table  2). 

Bed  armoring 

Bed  armoring  was  examined  in  two  model  runs  that  illustrated  alternate 
mechanisms  within  the  model.  In  the  first  run,  critical  shear  stress  for 
erosion  was  varied  with  depth  into  the  bed.  Critical  stress  was  specified  as 
2  dyne  cnr2  in  the  upper  two  layers  of  a  five-layer  bed.  Critical  stress  was 
12  dyne  cnr2  in  the  remaining  layers.  A  constant  shear  stress  of  10  dyne  cm-2 
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Figure  6.  Results  from  test  of  settling  through  the  water  column.  At  the  end  of  a  30-day  run, 
computed  concentrations  of  fixed  solids  (original  constituent),  LPOC  (modified  constituent), 
and  clay/silt  (new  constituent)  are  identical  in  three  layers. 


Table  2.  Transport  test  results.  Concentrations  are  presented  for  each  cell  in  the  30-box  test  grid  (Figure  1)  after 

a  30-day  model  run. 


Fixed  Solids  (g  nrv3) 

6.1252 

5.5968 

5.5136 

5.4992 

5.4958 

5.4882 

5.4442 

5.2479 

4.6358 

3.3559 

10.155 

9.6632 

9.5972 

9.5871 

9.5826 

9.5582 

9.4182 

8.8444 

7.2296 

4.2711 

17.3246 

16.7359 

16.6663 

16.6566 

16.6492 

16.5994 

16.3146 

15.167 

12.0135 

6.4678 

LPOC  (g  nr3) 

6.1252 

5.5968 

5.5136 

5.4992 

5.4958 

5.4882 

5.4442 

5.2479 

4.6358 

3.3559 

10.155 

9.6632 

9.5972 

9.5871 

9.5826 

9.5582 

9.4182 

8.8444 

7.2296 

4.2711 

17.3246 

16.7359 

16.6663 

16.6566 

16.6492 

16.5994 

16.3146 

15.167 

12.0135 

6.4678 

Clay/Silt  (g  nr3) 

6.12 

5.6 

5.51 

5.5 

5.49 

5.49 

5.44 

5.25 

4.63 

3.35 

10.15 

9.66 

9.6 

9.59 

9.58 

9.56 

9.42 

8.84 

7.23 

4.27 

1733 

16.74 

16.67 

16.66 

16.65 

16.6 

16.32 

15.17 

12.02 

6.47 
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was  imposed.  The  bed  was  initiated  with  20%  of  each  particulate  compo¬ 
nent  and  the  water  column  was  initiated  with  zero  suspended  particles. 
Two-dimensional  circulation  was  imposed  and  the  model  was  run  for 
30  days. 

The  upper  layer  (5  cm)  of  the  bed  eroded  into  the  water  column  imme¬ 
diately  (Figure  7).  Although  the  gross  erosion  rate  in  the  second  layer  was 
identical  to  the  first  layer,  a  portion  of  the  previously  eroded  material 
simultaneously  settled  so  that  net  erosion  (erosion  -  deposition)  slowed  and 
the  second  layer  took  much  longer  than  the  first  to  disappear.  Erosion 
ceased  after  the  upper  two  layers  disappeared  and  the  imposed  stress  no 
longer  exceeded  the  critical  shear  stress  for  erosion;  the  bed  armored.  An 
unexpected  result  was  the  subsequent  accumulation  of  material  in  Layer  3. 
Material  settling  into  Layer  3,  now  at  the  bed  surface,  acquired  the  proper¬ 
ties  of  Layer  3,  notably  critical  shear  stress.  The  thickness  of  this  layer 
increased  from  5  to  nearly  7  cm  as  material  accumulated;  Layer  2  was  not 
restored  to  replace  the  layer  that  previously  disappeared. 


Figure  7.  Results  from  test  of  bed  armoring.  Bed  Layer  3  armors  because  the  critical  shear 
stress  for  erosion  exceeds  the  imposed  shear  stress. 
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In  a  subsequent  test,  initial  conditions  in  the  bed  were  varied.  The  upper¬ 
most  layer  was  initialized  at  40%  clay/silt,  40%  fine  sand,  and  20% 
medium  sand  (no  LPOC  or  RPOC).  The  second  layer  was  initialized  at  10% 
clay/ silt,  10%  fine  sand,  and  80%  medium  sand.  All  layers  below  the  upper 
two  were  split  evenly  into  five  components.  The  vertical  distribution  of 
critical  shear  stress  and  the  imposed  shear  stress  were  retained  from  the 
first  test.  As  with  the  first  test,  the  upper-most  sediment  layer  rapidly 
eroded  into  the  water  column  (Figure  8).  Erosion  of  the  second  layer, 
however,  was  significantly  slower  than  in  the  first  test.  In  fact,  previously 
eroded  material  accumulated  in  Layer  2.  Subsequent  erosion  of  Layer  2 
approached  zero  despite  the  excess  of  imposed  stress,  10  dyne  cm  2,  over 
critical  shear  stress,  2  dyne  cm  2.  In  this  case,  the  bed  armored  because  it 
consisted  primarily  of  coarse  heavy  material  with  a  low  probability  of 
suspension.  The  initially  small  fractions  of  fine  material  decreased, 
causing  the  fraction  of  medium  sand  to  increase  from  80%,  at  the 
initiation  of  the  run,  to  88%  at  the  end.  Both  of  these  tests  indicated  the 
SEDZLJ  code  provided  acceptable  representation  of  bed  armoring. 


Figure  8.  Results  from  test  of  bed  armoring.  Although  imposed  shear  stress  exceeds  critical 
shear  stress,  bed  Layer  2  does  not  erode.  The  probability  of  suspension  for  medium  sand  is 

low. 


ER  DC/EL  TR-12-17 


14 


Continuous  erosion 

For  this  test,  a  five-layer  bed  was  set  up.  The  upper  four  layers  were  5  cm 
thick  while  the  lowest  layer  was  80  cm.  Critical  shear  stress  for  erosion 
was  specified  as  2  (upper  2  layers)  to  5  dyne  cm  2,  while  continuous  shear 
stress  of  20  dyne  cm-2  was  imposed.  The  upper  four  layers  eroded  into  the 
water  column  immediately.  The  bottom  layer  demonstrated  continuous, 
smooth  erosion  for  the  duration  of  the  30-day  run  (Figure  9).  This  test 
demonstrated  the  ability  of  SEDZLJ  to  respond  to  erosion  events  by  losing 
layers  and  by  eroding  continuously  from  a  single  layer. 


Figure  9.  Results  from  test  of  continuous  erosion.  The  upper  four  bed  layers  erode 
immediately.  The  bottom  layer  erodes  continuously  for  30  days. 


Continuous  deposition 

For  this  test,  shear  stresses  were  specified  so  that  erosion  would  not  occur. 
Initial  conditions  and  boundary  conditions  in  the  water  column  were 
specified  to  provide  a  large  reservoir  of  material  capable  of  settling.  Two- 
dimensional  circulation  was  enabled  and  the  model  was  run  for  30  days. 
Material  accumulated  continuously  in  the  surface  sediment  layer  (Fig¬ 
ure  10).  By  the  end  of  the  run,  Layer  1  had  increased  in  thickness  from  5  cm 
to  189  cm.  This  run  demonstrated  the  ability  of  SEDZLJ  to  continuously 
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accumulate  material.  This  run  also  reinforced  the  findings  from  the  tests  of 
bed  armoring;  the  thickness  of  the  surface  layer  increased  as  material 
accumulated  but  no  new  layer  was  created  despite  the  large  amount  of 
material  deposited. 


Figure  10.  Results  from  test  of  continuous  deposition.  The  surface  bed  layer  accumulates 
material  continuously.  No  new  layer  is  added  to  accommodate  the  additional  material. 
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4  Application  to  Prototype  System 

Lake  George,  Florida 

Additional  testing  of  the  code  was  gained  through  application  of  the 
combined  ICM/SEDZLJ  models  to  Lake  George,  Florida.  This  system  was 
selected  because  it  was  the  subject  of  a  previous  sediment  transport 
application  based  on  a  different  formulation.  Application  to  Lake  George 
facilitated  testing  of  SEDZLJ,  since  forcing  functions  and  several  necessary 
parameter  values  were  available.  The  previous  application  also  provided  a 
comparison  for  evaluating  the  capabilities  and  advantages  of  SEDZLJ. 

Lake  George  is  a  shallow,  freshwater  lake  (area  *  210  km2,  average  depth  « 
3.1  m)  situated  in  northeast  Florida  at  the  head  of  the  lower  St.  Johns  River 
(Figure  11).  Water  leaves  the  lake  through  a  forked  passage  at  the  northern 
end.  Flow  is  primarily  northwards,  towards  the  St.  Johns,  but  actions  of 
wind  and  tide  occasionally  force  water  back  from  the  St.  Johns  into  Lake 
George.  Lake  George  can  be  divided  into  three  regions,  based  on  computed 
currents  and  bottom  stress.  The  first  region  is  the  long,  narrow  neck  at  the 
northern  end  that  connects  the  lake  with  the  lower  St.  Johns  River.  In  this 
region,  currents  are  primarily  tidal  and  have  a  magnitude  of  10  cm  sA  The 
second  region  is  the  open  expanse  of  the  lake  itself.  Currents  in  this  section 
are  also  primarily  tidal  but  an  order  of  magnitude  less  than  in  the  northern 
neck.  The  third  region  is  the  tail  at  the  lower  end  that  connects  the  lake  to 
the  upper  St.  Johns  River.  In  this  tail,  the  tidal  signal  is  absent,  currents  are 
primarily  riverine,  and  vary  widely  in  magnitude.  Velocity-generated 
bottom  shear  stresses,  which  are  proportional  to  velocity  squared,  reflect 
and  exaggerate  the  current  regimes.  In  the  northern  neck,  shear  stresses  are 
periodic  and  have  a  magnitude  of  1  dyne  cm-2  (Figure  12).  Velocity¬ 
generated  stresses  are  periodic  in  the  open  lake  as  well,  but  are  one  to  two 
orders  of  magnitude  less  than  in  the  northern  neck  (Figure  13).  In  the  tail, 
velocity-generated  stresses  vary  over  orders  of  magnitude  and  are,  during 
flow  events,  the  largest  in  the  system  (Figure  14).  In  the  open  lake,  bottom 
stress  is  exerted  primarily  by  wind-generated  waves  (Figure  15)  while  in  the 
narrow,  convoluted,  neck  and  tail,  fetch  is  limited  such  that  wind-generated 
stress  on  the  bottom  is  zero.  At  their  largest,  wind-generated  stresses  in  the 
lake  are  much  less  than  the  velocity-generated  stresses  in  the  inlet  and 
outlet.  While  algorithms  exist  for  combining  shear  stresses  generated  by 
currents  and  wind,  the  disparity  in  the  magnitudes  of  these  forcing 
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Figure  11.  Lake  George,  Florida. 
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Figure  12.  Computed  velocity-generated  bottom  shear  stress  in  the  northern  neck  of 

Lake  George. 
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Figure  13.  Computed  velocity-generated  bottom  shear  stress  in  the  open  lake. 


Figure  14.  Computed  velocity-generated  bottom  shear  stress  in  the  southern  “tail”  of 

Lake  George. 
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Figure  15.  Observed  wind  and  computed  bottom  shear  stress  from  two  sources,  wind 
waves  and  currents,  in  the  open  lake,  January  1997. 

functions  in  the  different  regions  of  the  lake  leads  to  a  simplifying 
assumption:  shear  stress  in  the  neck  and  tail  is  derived  exclusively  from 
currents,  while  shear  stress  in  the  open  lake  is  derived  exclusively  from 
waves. 

The  Lake  George  model  operates  on  a  three-dimensional  grid  consisting  of 
563  cells  in  the  horizontal  plane  by  two  cells  deep  (Figure  16).  CE-QUAL- 
ICM  treats  each  cell  as  a  control  volume  and  applies  the  three-dimensional 
mass-conservation  equation  with  appropriate  sources  and  sinks.  The 
application  period  is  three  years,  commencing  on  January  1, 1997. 

Previous  Lake  George  application 

The  previous  Lake  George  application  incorporated  a  complete  eutrophica¬ 
tion  model  (Cerco  and  Noel  2009)  including  algal  production,  nutrient 
cycling,  and  sediment  diagenesis,  as  well  as  suspended  solids.  A  single  fixed 
solids  class  was  considered.  Computed  particulate  organic  carbon  was 
converted  to  volatile  solids  (VSS  =  2.5  *  POC)  and  combined  with  fixed 
solids  for  comparison  with  observed  total  suspended  solids  (TSS).  Algo¬ 
rithms  for  computation  of  fixed  solids  followed  Luettich  et  al.  (1990),  who 
modeled  suspended  sediments  in  a  shallow  lake  perturbed  by  wind  events. 
Settling  and  erosion,  the  sources  and  sinks  for  inorganic  suspended  solids, 
are  represented  by: 
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Figure  16.  Computational  grid  superimposed  on  an  aerial  photo  of  Lake  George  and 
surroundings.  The  lengthy  passage  (“neck”  at  top)  to  the  lower  St.  Johns  River  is 
included  in  the  grid  as  is  a  small  portion  of  the  major  lake  inflow  (“tail”  at  bottom). 


8C  S[C—Cb]  E 

——  =  -  ws — L— - H - 

8t  8z  A z 


(10) 


where 

C  =  inorganic  suspended  solids  concentration  (g  m  3) 

Cb  =  background  concentration  (g  nrs) 
ws  =  settling  velocity  (m  d1) 

E  =  erosion  rate  (g  nr2  d1) 

Az  =  thickness  of  grid  cell  adjoining  bottom 
z  =  vertical  coordinate  (m) 
t  =  time  (d) 

The  suspended  solids  model  incorporates  a  background  concentration  of 
fine  material,  which  does  not  settle. 

The  erosion  term  is  applied  only  in  cells  adjoining  the  bottom  (z  =  d)  and 
is  proportional  to  excess  shear  stress: 
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T —  TC 

E  =  ivs-Cref • -  (11) 

t ref 

where 

Cref  =  reference  concentration  (g  m  3) 

zc  =  critical  shear  stress  for  erosion  (dyne  cm-2) 
rref  =  reference  shear  stress  (dyne  cm  2) 
n  =  exponent  of  order  l 

The  model  is  grounded  on  the  assumption  of  long-term  equilibrium 
between  sediment  deposition  and  resuspension;  the  bed  solids  concentra¬ 
tion  is  effectively  constant.  Settling  velocity  is  explicitly  incorporated  into 
the  erosion  term.  This  relationship  necessarily  accompanies  the  assumption 
of  a  constant  bed  mass.  High  settling  velocity  must  be  accompanied  by  high 
erosion  rate,  and  vice-versa,  or  the  bed  would  not  be  in  equilibrium.  The 
reference  concentration  is  related  to  an  equilibrium  concentration  that 
would  be  attained  under  conditions  of  constant  excess  (x  >  xc)  shear  stress: 

n 

Ceq  =  Cref  -  — — ^  (12) 

r ref 

where 

Ceq  =  equilibrium  suspended  solids  concentration  (g  nrs) 

In  this  model,  solids  continuously  settle  from  the  water  into  the  bed;  there 
is  no  critical  shear  stress  for  deposition. 

Resuspension  of  organic  carbon 

As  originally  formulated,  the  representation  of  the  sediment  bed  in  the 
diagenetic  model  neglected  resuspension.  The  mass  balance  (DiToro  2001) 
included  net  deposition  from  the  water  column,  diagenesis  (decay),  and 
burial  to  deep,  inactive,  sediments: 


H~—  —  j_]c.HCs—wbCs 
dt 


(13) 
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in  which 

Cs  =  concentration  in  sediment  bed  (g  m  3) 

H  =  depth  of  sediment  bed  (m) 

J  =  net  deposition  rate  (=  wnet  •  C,  g  nr2  d1) 
k  =  diagenesis  rate  (d1) 
wb  =  burial  rate  (m  d1) 


For  Lake  George,  erosion  was  incorporated  into  the  diagenetic  sediment 
bed  by  defining  a  resuspension  velocity: 


R 


Es 

Ms 


(14) 


where 

R  =  resuspension  velocity  (m  d1) 

Es  =  inorganic  solids  erosion  rate  (g  nr2  d1) 

Ms  =  inorganic  solids  concentration  in  bed  (g  nrs) 

One  interpretation  of  the  resuspension  velocity  is  the  rate  of  change  of  the 
bed  thickness  if  the  bed  solids  concentration  remains  constant.  The  mass 
balance  for  the  bed  then  becomes: 

d 

H——  =  J—RCs—kHCs—wbCs  (15) 

dt 

and  the  erosional  source  to  the  water  column  is 

E  =  RCs  (16) 

Deposition  is  now  represented  by  the  settling  through  the  water  column, 
ws,  rather  than  the  previous  net  settling  rate. 

Parameter  set 

Parameters  for  the  sediment  resuspension  algorithms  were  empirically 
derived.  Initial  values  were  input  to  the  model  and  adjusted  in  subsequent 
model  runs  until  satisfactory  agreement  was  reached  with  observed  solids 
data  in  the  water  column.  This  practice  was  consistent  with  parameter 
evaluation  conducted  in  similar  studies  (Luettich  et  al.  1990;  Bailey  and 
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Hamilton  1997;  Hawley  and  Lesht  1992;  James  et  al.  2004a,  2004b)  and 
resulted  in  similar  values  (Table  3).  Two  key  parameters  in  the  resuspension 
algorithm  were  the  critical  shear  stress  for  erosion  xc  and  the  erosion 
potential  (equivalent  to  the  product  of  ws  and  Cref)  E.  The  critical  shear 
stress  for  the  neck  and  tail  was  three  times  the  value  derived  for  the  open 
lake,  indicating  the  bottom  sediments  of  the  neck  and  tail  are  more  resistant 
to  erosion  than  the  lake.  The  erosion  potential  of  the  lake  was  an  order  of 
magnitude  higher  than  the  neck  and  tail,  representing  a  reservoir  of  fine, 
easily  eroded  material  that  is  absent  in  the  high-velocity  regions. 


Table  3.  Parameter  set  from  previous  model  application. 


Lake 

Location 

ws, 

Cb, 

Cref, 

E, 

xc,  dyne 
crtT2 

-cref, 

dyne 

errf2 

n 

Reference 

met1 

gnr3 

gm-3 

g  nr2  d1 

Lake  George, 
FL 

Open  lake 

40 

2 

100 

4000 

0.5 

1 

1 

Lake  George, 
FL 

Neck  and  tail 

40 

2 

10 

400 

1.5 

1 

1 

Lake 

Balaton, 

Flungary 

One  of  three 
acceptable 
parameter  sets 

19 

15 

23.5 

447 

0.95 

0.88 

Luettich  et  al. 
1990 

Thomsons 

Lake, 

Australia 

One  of  two 
alternate 
parameter  sets 

2.6 

2 

766 

1992 

0 

1 

1 

Bailey  and 
Flamilton 

1997 

Lake  St. 

Clair,  Ml 

Range  of 
Parameters 
from  19  data 
sets 

1.7  to 
121 

1  to  4 

17  to 
158 

53  to 
3600 

0 

1 

1 

Flawley  and 
Lesht  1992 

Lake 

Christina, 

MN 

No-macrophyte 

conditions 

2.6 

100 

207 

538 

1.4 

0.67 

James  et  al. 
2004a 

Peoria  Lake, 

IL 

47 

120 

207 

9658 

2.3 

1 

0.67 

James  et  al. 
2004b 

SEDZLJ  application  to  Lake  George 

The  Lake  George  data  set  consists  primarily  of  TSS  and  POC  observations  in 
the  water  column,  along  with  the  composition  of  a  few  sediment  cores.  The 
previous  suspended  solids  model  was  selected,  in  part,  because  its  applica¬ 
tion  required  minimal  data  beyond  that  which  was  available.  Parameters  for 
which  observations  were  not  available  were  evaluated  through  recursive 
application  of  the  model  and  comparison  with  observed  properties.  SEDZLJ 
is  more  demanding  in  terms  of  required  parameters  and  application  should 
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be  accompanied  by  a  specific  data  collection  effort.  In  the  absence  of 
observations,  parameters  for  this  test  are  derived  from  the  previous  applica¬ 
tion,  when  possible,  or  default  values  provided  by  the  contractor  with  the 
code  are  employed.  No  optimization  of  the  SEDZLJ  parameters,  through 
recursive  model  application,  is  attempted.  The  selected  parameter  values 
suit  the  purposes  of  the  tests  -  to  evaluate  the  model  under  prototype 
conditions  and  to  compare  results  with  an  alternate  application. 

Settling  velocities 

Settling  velocities  for  inorganic  particles  are  computed  within  the  SEDZLJ 
code  based  on  specified  particle  diameters.  Diameters  for  the  three 
inorganic  classes  were  adapted  from  default  values  supplied  by  the 
contractor  (Table  4).  Settling  velocities  for  POC  were  characteristic  of 
values  used  in  eutrophication  models  (Table  4). 


Table  4.  Particle  diameters  and  settling  velocities  for 
SEDZLJ  application. 


Particle 

Diameter  pm 

Settling  Velocity  m  d1 

Clay/Silt 

5 

1.3 

Fine  Sand 

90 

396 

Medium  Sand 

500 

5069 

LPOC 

0.1 

RPOC 

0.1 

Bed  composition 

Analysis  of  the  sediment  cores  suggested  that  bed  sediments  in  the  neck, 
which  are  subject  to  regular,  high,  shear  stress,  were  denser  and  contained 
less  carbon  than  sediments  in  the  open  lake,  subject  to  irregular,  lower 
shear  stress.  The  sediment  bed  was  divided  into  five  layers.  The  upper  four 
layers  were  initially  5  cm  thick,  the  bottom  layer  was  10  cm.  Initial 
composition  of  the  upper  four  layers  was  based  on  observations  (Table  5). 
The  bottom  layer  was  assigned  composition  of  medium  sand  to  provide  a 
stable  base  resistant  to  erosion  in  the  event  the  upper  four  layers  eroded 
away.  Initial  bed  composition  in  the  tail,  also  subject  to  high  shear  stress, 
was  assigned  the  same  values  as  the  neck. 
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Table  5.  Initial  bed  characteristics  for  SEDZU  application. 


Characteristic 

Neck,  Tail 

Lake 

Bulk  Density,  upper  4  layers,  g  cnr3 

1.2 

0.4 

Clay/Silt  Fraction,  upper  4  layers 

0 

0.31 

Fine  Sand  Fraction,  upper  4  layers 

0.496 

0.303 

Medium  Sand  Fraction,  upper  4  layers 

0.496 

0.303 

LPOC  Fraction,  upper  4  layers 

0.004 

0.042 

RPOC  Fraction,  upper  4  layers 

0.004 

0.042 

Bulk  Density,  bottom  layer,  g  cm'3 

1.2 

1.2 

Clay/Silt  Fraction,  bottom  layer 

0 

0 

Fine  Sand  Fraction,  bottom  layer 

0 

0 

Medium  Sand  Fraction,  bottom  layer 

1 

1 

LPOC  Fraction,  bottom  layer 

0 

0 

RPOC  Fraction,  bottom  layer 

0 

0 

Critical  shear  stress 

Much  of  the  parameter  assignment  for  Lake  George  is  based  on  the 
division  into  regions  dominated  by  current -generated  shear  stress  (neck 
and  tail)  versus  regions  dominated  by  wave-generated  shear  stress  (open 
lake).  Critical  shear  stresses  for  surficial  sediments  in  these  regions  were 
adapted  from  the  previous  application  (Table  6).  Critical  shear  stress  was 
assumed  to  increase  with  depth  into  the  bed.  An  artificial,  high  value  was 
assigned  to  the  bottom  layer  to  prevent  erosion.  Critical  shear  stress  for 
deposition,  0.5  dyne  cm  2,  was  selected  from  the  range  of  values 
recommended  in  the  SEDZU  users’  guide. 


Table  6.  Critical  shear  stress  for  erosion,  SEDZU 
application. 


Layer 

Neck,  Tail,  dyne  cm  2 

Lake,  dyne  cnr2 

1 

1.5 

5 

2 

1.65 

0.55 

3 

1.82 

0.61 

4 

2 

0.67 

5 

999 

999 
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Erosion  rate 

SEDZLJ  calculation  of  erosion  rate  (cm  s1)  relies  on  two  parameters,  A  and 
n  (Equation  8).  Previously  determined  potential  erosion  rates  for  the 
system,  in  units  consistent  with  the  previous  model,  were  400  g  nr2  d1  for 
the  neck  and  tail  and  4000  g  nr2  d1  for  the  open  lake.  Division  of  these  rates 
by  bulk  density  (1.2  g  cnr3  in  the  neck  and  0.4  g  cm  3  in  the  lake)  followed  by 
units  conversion  yields  characteristic  values  for  parameter  A:  3.86  x  icr7  in 
the  neck  and  tail;  1.16  x  io_5  in  the  lake.  Parameter  n  was  assigned  a  value 
of  1,  consistent  with  the  previous  application.  The  exponent  in  the  cohesive 
sediment  hiding  factor  was  assigned  the  value  0.85  recommended  in  the 
SEDZLJ  users’  guide. 

Results 

The  SEDZLJ  application  was  “spun  up”  from  initial  conditions  for  a  period 
of  three  years.  Conditions  in  the  water  column  and  bed  at  the  end  of  the 
spin-up  were  used  as  initial  conditions  for  a  three-year  application.  The 
integration  time-step  was  1  minute. 

Time  series  comparisons 

Solids  observations  in  the  water  column  are  available  at  four  locations 
(Figure  17): 

•  MSJFGF  -  in  the  northern  neck,  subject  to  periodic  current-generated 
shear  stresses 

•  LG12  -  At  the  juncture  between  the  neck  and  open  lake,  shielded  from 
wind-generated  waves 

•  MSJLGM  -  At  the  center  of  the  open  lake,  subject  to  shear  stresses 
from  wind-generated  waves 

•  MSJLGS  -  At  the  southern  end  of  the  lake,  influenced  by  inflows  from 
the  riverine  tail 

TSS  comparisons  at  the  center  of  the  lake  reveal  characteristic  differences 
between  the  two  applications:  computed  TSS  values  in  the  previous  applica¬ 
tion  are  “spikey”  with  high  concentrations  that  settle  out  rapidly  while 
computed  TSS  in  SEDZLJ  attain  lower  concentrations  but  remain  in  the 
water  column  longer  (Figures  18  -  20).  TSS  computations  from  the  previous 
model  also  demonstrate  the  presence  of  the  2  gm  nr3  background  concen¬ 
tration,  which  is  absent  from  SEDZLJ.  The  longer  residence  time  of 
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particles  in  the  water  represented  by  SEDZLJ  can  be  attributed  to  the  lower 
settling  velocities  of  SEDZLJ  clay/silt  versus  fixed  solids  in  the  previous 
model  (1.3  m  d1  vs.  40  m  d1)  and  of  SEDZLJ  LPOC  and  RPOC  versus 
identical  particles  in  the  previous  model  (0.1  m  d1  vs.  40  m  d1).  Although 
SEDZLJ  sand  components  possess  higher  settling  velocities  than  fixed 
solids  in  the  previous  model,  computed  sand  fractions  are  negligible 
compared  to  clay/ silt  (Figure  21).  The  highly  variable  nature  of  the 
computations  and  the  scarcity  of  observations  at  this  location  prevent  a 
visual  judgment  of  the  accuracy  of  one  model  versus  the  other. 


Figure  17.  Lake  George  sample  stations. 
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Figure  18.  Computed  and  observed  TSS  at  station  MSJLGM  for  two  suspended  solids  models, 
SEDZU  and  an  “old”  model  based  on  long-term  equilibrium  between  deposition  and 

resuspension. 


Figure  19.  Computed  fixed  solids  at  station  MSJLGM  for  two  suspended  solids  models, 
SEDZU  and  an  “old”  model  based  on  long-term  equilibrium  between  deposition  and 

resuspension. 
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Figure  20.  Computed  and  observed  POC  at  station  MSJLGM  for  two  suspended  solids  models, 
SEDZU  and  an  “old”  model  based  on  long-term  equilibrium  between  deposition  and 

resuspension. 


In  the  northern  neck,  the  higher  concentrations  and  rapid  resuspension/ 
deposition  provided  by  the  previous  model  apparently  result  in  better 
representation  of  the  observations  than  from  SEDZU  (Figure  22).  It  should 
be  noted  that  the  previous  model  was  calibrated  to  match  the  observations 
while  no  adjustment  was  performed  to  the  initial  SEDZU  parameter  set. 
Presumably,  SEDZU  could  be  adjusted  to  improve  performance  in  this 
region.  The  characteristic  differences  between  the  two  representations 
noted  at  MSJLGM  prevail  at  LG12  (Figure  23)  and  MSJLGS  (Figure  24)  as 
well. 

Solids  budgets 

Although  there  are  solids  inputs  to  Lake  George  from  springs  and,  during 
storms,  from  distributed  flows,  the  preponderance  of  solids  enter  at  the 
southern  end,  from  the  Upper  St.  Johns  River.  The  only  exit  is  at  the 
northern  end,  into  the  Lower  St.  Johns.  For  practical  purposes,  lake  solids 
budgets  can  be  constructed  using  mass  fluxes  at  these  two  locations.  The 
budgeting  algorithms  in  ICM  were  used  to  create  monthly  (Figures  25  and 
26)  and  three-year  average  budgets  for  fixed  solids  and  POC  (Table  7).  For 
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both  model  versions,  fixed  solids  exports  exceed  fixed  solids  loads.  For 
SEDZLJ,  fixed  solids  outflow  exceeds  inflow  by  a  factor  of  six,  indicating 
net  erosion  from  the  sediment  bed.  The  excess  is  largely  in  the  form  of  fine 
sand;  a  small  import  of  medium  sand  is  computed.  For  the  previous 
model,  fixed  solids  outflow  exceeds  inflow  by  a  factor  of  30.  This  excess  is 
paradoxical,  since  the  model  is  constructed  around  the  concept  of  long¬ 
term  equilibrium  in  the  sediment  bed.  The  inclusion  of  settling  velocity  in 
the  erosion  term  is  intended  to  promote  bed  equilibrium. 


Clay,  MSJLGM  Fine  Sand,  MSJLGM 


Medium  Sand,  MSJLGM 


POC,  MSJLGM 


Figure  21.  Computed  individual  sediment  components  from  SEDZLJ  at  station 

MSJLGM. 
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TSS,  MSJFGF 


- sedzlj 

0  observed 

100 


0  365  730  1095 

days  since  Jan  1,  1997 


Figure  22.  Computed  and  observed  TSS  at  station  MSJFGF  for  two  suspended  solids  models, 
SEDZU  and  an  “old”  model  based  on  long-term  equilibrium  between  deposition  and 

resuspension. 


Figure  23.  Computed  and  observed  TSS  at  station  LG12  for  two  suspended  solids  models, 
SEDZU  and  an  “old”  model  based  on  long-term  equilibrium  between  deposition  and 

resuspension. 
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Figure  24.  Computed  and  observed  TSS  at  station  MSJLGS  for  two  suspended  solids  models, 
SEDZU  and  an  “old”  model  based  on  long-term  equilibrium  between  deposition  and 

resuspension. 


Figure  25.  Monthly  fixed  solids  budget  (kg  d1)  from  two  models. 
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Table  7.  Solids  budgets  for  two  models,  averaged  over  the  simulation  period. 


Budget 

SEDZLJ 

Previous  Model 

Inflow 

Outflow 

Inflow 

Outflow 

Clay/Silt,  kg  d-1 

8,794 

20,923 

Fine  Sand,  kg  d  1 

8,798 

132,336 

Medium  Sand,  kg  d 1 

9,152 

-714 

Fixed  Solids,  kg  d1 

24,208 

152,546 

24,208 

786,652 

POC,  kg  d1 

10,390 

45,882 

16,175 

24,230 

POC  outflow  computed  by  SEDZLJ  exceeds  inflow  by  roughly  the  same 
factor  as  fixed  solids,  once  again  indicating  net  erosion  of  the  bed.  The 
previous  model  POC  inflow  and  outflow  are  roughly  in  balance,  in  contrast 
to  fixed  solids.  The  contrast  is  attributed  to  modeling  of  the  complete 
carbon  cycle,  which  includes  mass  balance  in  the  bed.  These  results 
indicate  that  rigorous  sediment  transport  modeling  must  include  bed  mass 
balance.  Although  the  present  SEDZLJ  representation  of  Lake  George  is 
out  of  balance,  the  formulation  provides  for  balance  to  be  achieved 
through  additional  tuning  of  the  model,  ideally  aided  by  data  collection. 
The  previous  model  does  not  provide  this  opportunity  although,  perhaps, 
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whole-lake  mass  balances  could  be  employed  in  concert  with  parameter 
evaluation  to  obtain  a  consistent  model  using  the  simpler  formulation. 

Negative  sand 

Computation  of  negative  sand  concentration  was  noted,  in  particular  at  365 
days  into  the  simulation  (Figure  21).  Near  this  time,  quantities  of  sand  were 
eroded  from  the  riverine  portion  by  a  flood  event  and  deposited  in  the  lake. 
Apparently  numerical  instability  caused  negative  concentrations.  The 
settling  velocity  of  medium  sand  is  5069  m  d1  or  5.9  cm  s1.  Medium  sand 
settles  through  a  water  column  layer  of  average  thickness  in  25  s  while  the 
integration  time-step  is  60  s.  The  rapid  settling  produced  negative  concen¬ 
trations,  which  persisted  after  the  event  passed.  Experiments  were 
attempted  in  which  the  mass  allowed  to  settle  in  a  time-step  was  restricted 
to  the  mass  in  a  layer,  but  negative  concentrations  persisted.  In  reality,  this 
model  is  not  intended  to  transport  medium  to  coarse  sand.  The  role  of 
medium  sand  in  the  simulation  is  to  provide  a  substrate  of  stable  material  in 
the  bed.  Sand  should  be  transported  primarily  as  bedload,  which  is  absent 
from  this  model  version.  The  negative  sand  concentrations  are  harmless  but 
must  be  noted. 


Statistical  summary 

Statistical  summaries  were  prepared  for  TSS  and  POC  computed  by  both 
models  and  compared  to  observations.  In  keeping  with  summaries  of 
previous  ICM  simulations,  computed  statistics  included  mean  difference 
(MD): 


Y(P-O) 
=  - - 

n 


(17) 


where 

O  =  observation 

P  =  computation 

n  =  number  of  observations 


and  absolute  mean  difference  (AMD) 


P-0\ 


AMD 


n 


(18) 
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Positive  MD  indicates  the  model  computations  exceed  the  observations,  on 
average.  Negative  MD  indicates  the  model  computations  are  less  than 
observations,  on  average.  The  AMD  statistic  characterizes  the  magnitude 
of  the  difference  between  individual  computations  and  observations, 
regardless  of  sign.  Correlations  between  computations  and  observations 
were  calculated,  as  well. 

Statistical  summaries  indicate  the  two  models  are  equivalent  in  their 
representation  of  the  observations  (Table  8).  Both  models  under-predict 
TSS  by  the  same  amount,  «  6  g  m-3,  and  over-predict  POC  by  «  l  g  nr  3. 
Correlation  coefficients  between  computations  and  observations  are 
effectively  zero.  The  equivalence  of  the  two  computations  is  interesting  in 
view  of  the  two  different  representations  of  the  system.  The  previous 
model  described  a  system  with  a  constant  background  solids  concentration 
and  rapid  resuspension  and  deposition.  SEDZL  describes  a  system  with  no 
background  concentration  and  slower  resuspension  and  deposition.  Lack 
of  correlation  indicates  the  difficulty  in  exactly  matching  sparse,  sporadic, 
observations  with  the  model.  Rigorous  solids  computation  in  this  system 
requires  a  comprehensive  observational  program  including  field  samples 
and  sediment  characteristics.  Nevertheless,  the  goals  of  testing  SEDZLJ  in 
a  prototype  application  and  of  comparing  results  with  an  alternative  model 
have  been  fulfilled. 


Table  8.  Statistical  summary  of  computations  from  two  models. 


Component 

Previous  Model 

SEDZLJ 

TSS,  mean  difference,  g  nr3 

-5.2 

-6.9 

TSS,  absolute  mean  difference,  g  nr3 

9.8 

10 

POC,  mean  difference,  g  rrr3 

0.3 

1.6 

POC,  absolute  mean  difference,  g  nr3 

1.6 

2.7 
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5  Summary  and  Conclusions 

The  preceding  chapters  described  the  formulation  of  the  SEDZLJ 
sediment  transport  model,  coupling  of  the  model  with  the  ICM 
eutrophication  model,  validation  of  the  combined  codes  on  a  simplified 
grid,  and  application  of  the  combined  codes  to  a  prototype  system. 

Code  validation  tests  indicated  that  the  codes  are  correctly  coupled.  The 
combined  codes  met  specified  performance  criteria,  including: 

1.  Mass  conservation  in  water  and  bed. 

2.  Limited  sensitivity  to  variations  in  model  time-step. 

3.  Settling  of  new  ICM  state  variables  through  model  water  column  agrees 
with  settling  of  original  ICM  state  variables. 

4.  Transport  of  new  ICM  state  variables  agrees  with  transport  of  original  ICM 
state  variables. 

5.  Armoring  of  bed. 

6.  Smooth  erosion  of  bed. 

7.  Reasonable  behavior  of  the  model  for  accumulation  in  the  bed  with  no 
erosion. 

Additional  evaluation  and  testing  were  performed  through  application  of 
the  combined  models  to  a  prototype  system,  Lake  George,  Florida.  The 
code  performed  reasonably  through  a  three-year  simulation.  The  tests 
indicated,  however,  that  realistic  application  of  SEDZLJ  requires  an 
extensive  set  of  field  and  laboratory  observations. 

One  characteristic  revealed  by  the  prototype  simulations  was  the 
computation  of  negative  sand  in  the  water  column.  The  negative 
concentrations  apparently  result  from  numerical  aberrations  caused  by  the 
rapid  settling  of  sand  through  the  model  water  column.  Stability  criteria 
indicate  model  time-steps  «  10  s  might  be  required  to  properly  compute 
settling  of  sand.  In  reality,  the  present  model  is  not  intended  to  transport 
medium  to  coarse  sand,  since  bedload  is  not  represented.  The  role  of 
medium  sand  in  the  model  is  primarily  to  anchor  the  bed. 
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The  integration  time-step  in  the  present  model  was  reduced  to  60  s 
versus  15  min  for  application  of  the  ICM  eutrophication  algorithms  to  the 
prototype  system.  As  a  result  of  the  reduced  time-step  and  additional 
computational  demands  imposed  by  the  SEDZLJ  bed  model,  the  computa¬ 
tion  time  of  the  combined  ICM/SEDZLJ  models  increases  tremendously 
over  the  basic  eutrophication  model.  A  three-year  simulation  of  Lake 
George  eutrophication  consumes  4  hr  on  a  desktop  PC,  while  a  three-year 
simulation  of  suspended  solids,  using  SEDZLJ,  consumes  four  days.  This 
increase  is  despite  the  fact  that  the  eutrophication  simulation  includes 
26  state  variables,  while  the  suspended  solids  simulation  includes  only  six. 
Although  computation  time  can  be  reduced  by  moving  to  a  faster  computer, 
the  practical  application  of  the  combined  codes,  in  their  present  state,  is 
severely  limited.  Practicality  will  be  further  limited  when  computation  of 
toxicants  is  added  to  the  combined  codes.  Efforts  to  convert  the  combined 
codes  to  parallel  computation,  via  domain  decomposition,  are  presently 
underway.  Simultaneously,  ICM/SEDZLJ  is  being  combined  with  the 
RECOVERY  toxicant  model  (Boyer  et  al.  1994),  after  which  a  second  round 
of  parallelization  is  anticipated. 
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